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1. Introduction 

Since the sine, cosine, exponential integrals and hyperbolic sine, and cosine integrals are 
frequently encountered together in physical problems and their expansions have terms in common, 
we have incorporated these functions into Fail III. (For Parts I and II. sec. 1 ). 

While accuracy over the entire domain of definition remains our main concern, we have tended 
toward methods that also ensure efficiency, portability, and ease of programming and modification. 
The number of terms in series, the number of convergents in an iterative 4 process, the starting 
arguments for different methods, are all determined by the program as a function of word length, 
arguments, accuracy desired, etc. More realistic results are returned when error conditions are 
encountered. The proper analytic behavior of the function will always be retained to fin titer ensure 
correct limiting values, in particular of related functions and for purposes of differentiation and 
integration. 

In Parts I and II in addition to the implementing ANS FORTRAN program, we had included a 
driver (test) program and its results. Since either of these driver programs can be readily modified to 
compute other functions, we have omitted the driver program and in place of its results have 
included a table of correct results to 35 significant figures covering essentially the functional range of 
present computers. 



2. Mathematical Properties 

Relevant formulas are collected here for completeness and ease of reference. In keeping with 
the convention of the Handbook [l], 2 x here is a real variable. 



1 Automatic Computing Methods for Special Functions. Part I. Error, Probability, and Related Functions, I. Res. Nat. Bur. Stand. (U.S.), 74B. (Math. Sci), 
No. 3, 211-224 (July-Sept. 1970). Automatic Computing Method for Special Functions. Part II. The Exponential Integral E n (x). J. Kes. Nat. Bur. Stand. (U.S.), 
78B. (Math. Sci), No. 4, 1W-2H) (Oct.-Dec. I ( >7t). 

2 Figures in brackets indicate the literature references on page 
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A. Definitions 



sin t . 

Si(x) — at 

Jo t 



r*. x i f'cOSf - 1 . 

Li{x) = y + In * + or 

Jo * 



££(*) = - j — dt = \ -dt (x > 0) 

J _, * .L. t 

(For — dt = Ei(x), often denoted by -Ei(-x), see Part II.) 



sinh e £t(#) 4- E x (x) 

Sni{x) = at = 



„,., . , f cosh^-1 EiU) -£iU) 
Lni(x) = y + In x + at = 

Jo t 2 



y (Euler's constant) - 0.57721 56649 . . .. 



B. Series Expansions 

& - w = y ( - 1)A ^ + ' 

-o (2& + 1)(2* + D! 



Ci(x) = y + In x + > — ■ — - 

fc ti (2ftX2*)! 



£i(x) = y + In x + V — — - (x > 0) 

S (AX*)! 



shux) = y _ 

A (2A: + 1)(2* + 1)! 



Chi(\) = y + In x + > 



fc tl (2AX2A)! 



C. Continued Fraction 
-Ci(x) + i[S*(ac) - tt/2] = e"" [- — — : — — — ] (n < x) 

IX + 1 + IX + 1 + IX + . . . 

= E t (ix) 
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where f(x) 



D. Asymptotic Expansions 

Si{x) = 77/2 —f{x) cos x — g{x) sin x 
Ci{x) — f{x) sin x— gix) cos x 
1 * (-1)*(2A)! 



, , . 1 " (~l) fr (2A + 1)! 



— 2fr+l 

X £ = () A. 



€ J ' X ft 1 

£*(*) — ^ — (x > 0) 



Shi(x) = - [p(x) cosh x + g(x) sinh jc] 



Chi(x) = - [p(#) sinh .r + <?(#) cosh jt] 



v (2ft)! 

where p(x) - 21 — £7~ 

fr=0 * 



v (2A ~ + 1>! 

and 9 (x) - 1 , 

fc=0 X 



E. Special Values 

S£(0) = Ci(0) = - x £7(0) = - x 

Shi(0) = CAi(O) = - x 



F. Symmetry Relations 

(Jt > 0) 

Si{—x) = —Si(x) Ci(—x) = Ciix) — in 

Shi(-x) = —Shi(x) Chi(-x) = Chi(x) — in 



G. Interrelations 



Si(x) = — [E t (ix) - E^-ix)] + nil 
2 1 



Ciix) = - - [E^ix) + E x (—ix]\ 



Ei(x) = - - [E x (-x + £0) + EA-x - £0)] U > 0) 
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H. Value at Infinity 

lim Si{x) = 77/2 



I. Related Function 
Logarithmic Integral 



li(x) = Ei(h\(x)) 



(x > 1) 



3. Method 

Evaluation of the integrals by means of quadrature formulas suited to the particular type of 
integrand tends to be inefficient and inaccurate. For Si(x) and Ci(x), the use of the asymptotic 
expansion is not valid for moderate values of x, while the use of the continued fraction is inefficient 
and also inaccurate for small values of x. An examination of the series expansion for the functions 
indicates several difficulties. Summation of the alternating series for Si(x) and Ci(x) will lead to 
greater round-off errors as x increases. The partial sum at a particular value of k may be zero. 
Additionally there may be cancellations in adding the logarithmic term and/or Euler's constant for 
Ci(x), Ei(x) and Chi(x). The more rapidly accumulating round-off errors, in particular when 
summations are limited to a single register, eliminate the prolonged use of the series expansion. 
Since the maximum of Ci(x) occurs at ttI2 and Si(x) = tt/2 at x — 1.92, testing indicates x = 
2(PSLSC) as a reasonable starting point for the use of the continued fraction. The starting point for 
the valid use of the asymptotic expansion for Si(x) and Ci(x) does not coincide with the starting point 
for Ei{x) (Shi{x) and Chi(x)). Testing also indicates that fewer terms are needed in the continued 
fraction than in the asymptotic expansion. 

The following table gives an indication of the number of terms needed to obtain maximum 
machine accuracy for particular values of x with the various methods of computation. Throughout 
the paper, NBM is the maximum number of binary digits in the mantissa of a floating point number, 
and TOLER = 2~ NBM . 



Method 






Number 


f Terms 








NBM = 27 


NBM = 60 


x = 2 


Si(x) 


Ci(x) 


Ei(x) 


Si(x) 


Ci(x) 


Ei(x) 


Power Series 


7 


7 


13 


12 


12 


23 


Continued Fraetion 


24 


24 


— 


106 


106 


— 


(Even Form) 














Numerical Integration 


64 


128 


64 


— 


— 


— 


(Trapezoidal or Simpson's Rule) 














x = 24 














Power Series 


— 


— 


55 


— 


— 


80 


Asymptotic Expansion 


a 7,ll 


7,11 


13 


— 


— 


— 


Continued Fraction 


5 


5 


— 


14 


14 


— 


Numerical Integration 


512 


512 


1024 


— 


— 


— 


x = 48 














Power Series 


— 


— 


— 


— 


— 


119 


Asymptotic Expansion 


4,5 


4,5 


8 


16,23 


16,23 


31 


Continued Fraction 


4 


4 


— 


9 


9 


— 


x = 88 














Asymptotic Expansion 


3,4 


3,4 


6 


9.10 


9.10 


17 


Continued Fraction 


5 


5 


— 


8 


8 


— 



a We indicate the number of odd and even terms of the respective series. 



The most accurate, efficient, automatic methods for Si(x) and Ci(x) then are the power series 
and the continued fraction; for Ei(x), Shi(x) and Chi(x) the power series and the asymptotic 
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expansion. The lower limit (AELL) for the use of the asymptotic expansion may be shown to 
approximate |ln TOLER| = NBM(ln 2), where TOLER is the requested upper limit for the relative 

error. With this choice of the lower limit, one can also show that Shi(x) — Chi{x) — - Ei(x). It is 

necessary then to consider only the asymptotic expansion for Ei{x). 

The series computations have been so arranged that the maximum number of functions- may be 
obtained in a minimum of time. The even and odd terms of the series are summed independently 
both with and without the factor (— 1)*. Since Si(x) -tt/2 and Ci(x) are the imaginary and real parts 
respectively of the continued fraction expansion for Ei(ix\ then 4 would be a saving in computing 
time with options on the functions to be computed. Invalid results are initially supplied for all 
functions. With the parameter IC = I, Si(x) and Ci(x) only are computed; with IC = 2, Ei(x) and 
e~ x Ei(x) only; with IC = 3, Ei(x\ e J Ei(x\ Shi(x) and Chi{x) only and with IC — 4, all functions arc 
computed. 

The implementing program checks the input parameters. If IC is outside the range 1 — 1. the 
working indicator IND is automatically set equal to 4. Since Ei(x) is defined for positive 1 x only, if IC 
— 2 and x < 0, there is an error return and the indicator IERR is set equal to 1. II ' x < 0, IC = 3, 
Shi(x) and Chi(x) are computed; for IC = 4, Si(x) and CiU) are also computed, invalid results are 
returned for Ei(*) and e~ x Ei{x) and the indicator [ERR is set equal to 1. For x > 0, the indicator 
IERR is set equal to zero and valid results are returned only for the functions requested by the 
parameter IC (or IND). 

The computations are performed for positive x (=T) only. Before the return, use is then made 
of the symmetry relations. Various cases are treated independently. Among these are x = and x 
equal to or greater than the supplied upper limit (ULSC) for the sine and cosine routine. The 
appropriate values of the functions are supplied. To avoid unnecessary computations of the 
exponential function and possible overflows and underflows in the final results, if T is equal to or 
less than the upper limit for the relative error, the exponential of half the argument is set equal to 1. 
If T is equal to or greater than the computed limiting argument (XMAXHF) for Shi(x) and Chi(x), 
the maximum machine value (R1NE) is supplied for the exponential of half the argument as well as 
for Shi(x) and (Jhiix). The computed limiting argument (XMAXEI) for Eiix) can be shown to be 
approximately In RINE + ln[ln RINE + ln(ln RINF)] -1/ln RINK The value of XMAXHF is 
approximately XMAXEI + In 2. The computation of e xl2 provides a slight improvement in accuracy 
and an extension of the range of x. Throughout the program, overflows are avoided as well as 
underflows affecting accuracy. Other underflows are assumed to be set equal to zero. 

For |*|^ PSLSC ( = 2), all functions are computed by means of the power series. For 
\x\ > PSLSC, Si(x) and Ci{x) are computed by means of the continued fraction. Only if IC = 4 is 
the working indicator IND set equal to 3. The functions Ei(x), e~ x Ei(x), Shi(x) and Chi(x) are then 
computed by means of the power series or the asymptotic expansion depending on whether 
|*| < AELL or > AELL respectively. With NBM - 27, AELL - 18.7 and with NBM = 60, AELL 
~ 41.6. To avoid underflow, \x\ is tested against a lower limit argument PSLL(=2 V / AM IN). To 
simplify computation, AMIN, a minimum machine value is computed as the reciprocal of the 
maximum machine value (RINF). If \x | < PSLL, only the first term of the series of odd terms is 
used. 

The following series definitions are in use 

Si( x ) = SI = SUMS = V SGN (RKfTM(RK) IP = -I (RK = 1,3, . . .) 

Ci(x) = CI = SUMC 4- A LOG 4- EULER 

where SUMC = V SGN (RK)*TM(RK) IP = 1 (RK = 2,4, . . .) 

Ei(x) = EI = SUMET + SUMOT 4- X LOG 4- EULER 

Shi(x) = SHI = SUMOT - V TM(RK) IP = -1 (RK = 1,3, . . .) 

Chi(x) = CHI - SUMET[= v TM(RK)] + X LOG + EULER IP = 1 (RK = 2,4,. . .) 

with SGN(l) = 1, SGN(RK + 1) = -SGN(RK) for RK = 1,3,. . ., and SGN(RK 4- 1) = SGN(RK) 
for RK = 2,4, .... The term TM(RK) = [T k ikl]/k = PTM( RK)IRK where PTM(l) = T and 
PTM(RK 4- 1) - PTM(/xA)[77(RK + 1)], RK > 1. 
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The series of even and odd terms are always computed together. If the relative error RE 
computed as FM/|SUM| is less than the prescribed tolerance both series are considered to have 
converged. If IND = 1 or 4, SUM is replaced by SUMS or SUMC; otherwise by SUMET or 
SUMOT. To avoid underflow, in generating the terms for |x| < 2, if PTM < AMIN (RKfIT, the 
series are likewise considered to have converged. If the sum of terms is zero, the relative RE is 
automatically set equal to the maximum machine value. This condition is not encountered if the 
power series tor Si [x) and Ci{x) is restricted to the region \x\ < 2. It has been retained to permit the 
program's use for experimental purposes. 

To enable the continued fraction computations to be performed in double precision, since 
complex quantities are involved, the real notation only is used. Testing has also confirmed the 
improved accuracy and efficiency of this course. The continued fraction for Si{x) and Ci(x) in its 
"even" form 



E l (ix) = -Ci(x) + i[Si(x) - tt/2] = e~ u 



1 1 



>d(l 



\-\-ix — 3+ix — 5+ix — . . . 
e~ iJ [F] 

is evaluated in the forward direction. The first convergent F 1 /Gi = A X IE X where A\ = 1, A M = —{M 
- I) 2 , B M = 2M - 1 + ix. If we define 

F_, = 1, F = 0, G_i = and G = 1 

then successive convergents F M /G M for M = 1,2, ... may be obtained by the following recurrence 
relation 

Fm = BmFm-\ + A M F M - 2 

(*m = B m Gm-i + A M G m-2 

The continued fraction is considered to have converged either if in effect the relative error is equal 

to or less than the prescribed tolerance or the relative error increases. 

Since the successive convergents are complex, (RE) 2 is compared with (TOLER) 2 where (RE) 2 

F IC 

— — — —7 — ) 2 . Throughout the computation, to avoid overflow, there is scaling by the 

F M 'G m 

absolute maximum (=TMAX) of the real and imaginary parts of the numerator and denominator of 
the successive convergents. In addition, there is scaling by |TMAX| if the product of the real part of 
{B M — Am) and |TMAX| is equal to or greater than 1/4 the maximum machine value. 

The successive terms of the asymptotic expansion are likewise obtained by recurrence with T 
= 1 and T K = [A7jT]7#-i for A > 1. Since the sum of terms for Ei{x) is always greater than one, 
the term itself is a good approximation to the relative error. The summation is terminated when a 
term is less than the prescribed tolerance or the term is equal to or greater than the preceding term. 
In the latter case, the preceding term is subtracted from the summation to minimize the error. 



4. Range 

The range for Si(x) and Ci(x) (as well as the accuracy) is limited to the range (and accuracy) of 
the sine and cosine routine (\x\ < ULSC). For the UNIVAC 1108, namely, x < 2 21 in single 
precision and x<2 56 in double precision. For the function Ei(x), the range of x is essentially the 
range of the exponential routine. The function Ei{x) is set equal to the machine maximum (RINF) for 
x beyond XMAXEI, approximately 92.5 in single precision and 715.6 in double precision. For the 
function e~ x Ei{x) beyond x = ULSC only the first two terms of the asymptotic expansion are used. 
The functions Shi(x) and Chi(x) are set equal to the maximum machine value for x beyond 
XMAXHF, approximately 93.2 in single precision and 716.3 in double precision. 
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5. Accuracy and Precision 

Using the UNIVAC 1108 to compute the functions, the maximum relative error, except for 
regions in the immediate neighborhood of zeros, is 4.5 (-7) for single precision computations and 7.5 
(-17) for double precision computations. Various auxiliary functions are available to greater accuracy 
at intermediate points in the subroutine. For example, since Si (x) — > 77/2, Si{x) —77/2 should be 

taken as the imaginary part of the continued fraction. The functions Ci(x), Ei(x) and Chi(x) —y — In 
x are available from the sum of the appropriate series. 

The precision may be set lower than the maximum by varying the value of NBM or deleting NBM 
and setting a precomputed value of TOLER. The above relative errors give an indication of the 
allowance for round-off errors. 



6. Timing — UNIVAC 1108 Time/Sharing Executive System 

(The time estimates given below are highly dependent on the operating system environment and 
consequently should not be relied on for critical timing measurements.) 



Single Precision 

NBM = 27 



Double Precision 
NBM = 60 







For Si{x) and Ci{x) 






Region 


Time 
(seconds) 


Region 


Time 
(seconds) 


Method 


0(0.01)2 (201 values) 
2(.5)100 (197 values) 
Maximum Time/Evaluation 
(x = 2) 


0.40 
.56 
.0023 

.0093 


0(.01)2 (201 values) 
2(.5)10O (197 values) 

(x = 2) 


0.98 
2.06 
0.0059 
.070 


Power Series 
Continued Fraction 
Power Series 
Continued Fraction 



For Ei{x) 



0(.1)18 (181 values) 


0.54 


0(.2)41 (206 values) 


2.05 


Power Series 


18(.1)41 (231 values) 


.28 






Asymptotic Expansion 


41(. 25)100 (237 values) 


.25 


41(. 25)100 (237 values) 


0.70 


Asymptotic Expansion 


Maximum Time/Evaluation 


.0044 




.016 


Power Series 


(x = 18) 




(x = 41) 








.0015 




.0040 


Asymptotic Expansion 



7. Testing 

The double precision results obtained were compared against available published values. Check 
values were obtained, where appropriate, by overlapping the power series with either the asymptotic 
expansion or the continued fraction. Various forms of the continued fraction were also employed as 
well as numerical integration. Two multi-precision packages 3 were also utilized with varied precision. 
The single and double precision results agreed with the multi-precision results within the reported 
accuracy. 



8. Many-Place Tables 

In the appendix, we have included three tables; one for Si (x) and Ci(x), one for Ei(x) and 
e~ x Ei(x) and one for Shi(x) and Chi{x). The functions are tabulated to 35 significant figures for x = 
0, 10 J (10 7 ) 10- 7+1 with J = -2(1)2. 



•'(Private Communication) Peavy, Bradley A, A Multi-Precision Arithmetic Package for Use with the UNIVAC 1108. 
Wyatt, W. T. Jr., Lozier, D. W. and Orser, D. J., A Portable Extended-Precision Arithmetic Package and Library 
TOMS, Sept. 1975. 



ith FORTRAN Precompiler, ACM 
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9. Special Values 

Zeros 

Si(x s ) = W2 

x = 1.92644 7660 
*! = 4.89383 5953 
x 2 = 7.97268 2624 
Ci(x s ) = 

x = 0.61650 5486 
x x = 3.38418 0423 
* 2 = 6.42704 7744 
Ei(x) = 

x = 0.37250 74107 81366 63446 19918 66580 11913 
Chi(x) = 

x =0.52382 25713 89864 



Si(ir) = 1.85193 70519 82466 
Ci(tt/2) = 0.47200 06514 39569 



Si(2ir) = 1.41815 15761 32628 
Ci(3tt/2) = -0.19840 75606 92358 



Mc 



Minima 



Related Constants 



y (-1)" 

N % (2AT+1X2/V+1)! 



= S»(1) 



0.94608 30703 67183 01494 13533 13823 17965 



= Ci(l)-y 



N - (2N\2N)\ 



1 



ivti o 



y 

aS) (27V+1X27V + 1)! 



cm 

-Ei(l) -y 
Ei{\) 

--Shi(l) 



-0.23981 17420 00564 72594 38658 86193 25166 
.33740 39229 00968 13466 26462 03889 15076 

1.31790 21514 54403 89486 00088 44249 23183 
1.89511 78163 55936 75546 65209 34331 63426 

1.05725 08753 75728 51457 18423 54895 87795 



1 



*ti (2N)(2N)l -. 

y(Euler's constant) 

tt/2 

log P 2 



Chi(l)-y = 0.26065 12760 78675 38028 81664 89353 35387 

Chi(l) = .83786 69409 80208 24089 46785 79435 75630 

= 0.57721 56649 01532 86060 65120 90082 40243 

= 1.57079 63267 94896 61923 13216 91639 75144 

- 3.14159 26535 89793 23846 26433 83279 50288 

= 0.69314 71805 59945 30941 72321 21458 17656 
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Typical Tolerances and Their Natural Logarithms 

2-24 = 0.59604 64477 53906 25 (-7) 

2-27 = 74505 30595 9 23 8 2 8125 (-8) 

2" 36 = .14551 91522 83668 51806 64062 5 (-10) 

2" 48 = 35527 13678 80050 09293 55621 33789 0625 (-14) 

2~ 56 = .13877 78780 78144 56755 29539 58511 35253 90625 (-16) 

2- 60 = .86736 17379 88403 54720 59622 40695 95336 91406 25 (-18) 

2- 108 = .30814 8791 1 01957 73648 89564 70813 58837 09660 96263 

71446 21112 38390 20729 06494 14062 5 (-32) 

log e (2- 24 ) - -16.63553 23334 38687 42601 35709 14996 23763 

log e (2~ 27 ) = -18.71497 38751 18523 35426 52672 79370 76733 

log e (2- 36 ) = -24.95329 a5001 58031 13902 03563 72494 35645 

log 6 (2- 48 ) = -33.27106 46668 77374 85202 71418 29992 47526 

log e (2" 56 ) = -38.81624 21113 56937 32736 49988 01657 88781 

log e (2- 6() ) = -41.58883 08335 96718 56503 39272 87490 59408 

log e (2- 1()H ) = -74.85989 55004 74093 41706 10691 17483 06935 

Maximum and Minimum Machine Values and Their Natural Logarithms 
NBC = Number of binary digits in the (biased) characteristic of a floating point 
number 

2-<2 v/ " l +n < x < 22 v/ ' r -«-i 

NBC - 8 

2 127 = 0.17014 J 1831 60469 23173 16873 03715 88410 (39) 
2" 129 = . 14693 67938 52785 93849 60920 67152 78070 ( -38) 
loge(2 127 ) = 88.02969 1931 1 13054 29598 84794 25188 42414 
log,42 -129) = -89.41598 62922 32944 91482 29436 68101 77728 

NBC = 1 1 

2 102:* = 0.89884 65674 31 157 95386 46525 95394 51236 (308) 
2-1025 = 27813 42323 13400 17288 62790 89666 55050 (-308) 
log^2 1()2:i ) = 709.08956 57128 24051 53382 84602 51714 62914 
log e (2- 1025 ) = -710.47586 00739 43942 15266 29244 94630 98227 
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APPENDIX 

IMPLEMENTING PROGRAM 
LANGUAGE. AMERICAN NATIONAL STANDARD FORTRAN 
DEFINITIONS. Xr A REAL VARIABLE 

SKX) =INTEGRAL(SIN T/T)QT FROM TO X 
SI(-X)=-SI(X) 

CI(X) =GAMMA+LN X+INTEGRAL < (COS T-1)/T)DT FROM TO X 
CI(-X)=CI(X)-I PI 

EKX) =-P.V.INTEGRAL<EXP(-T)/T)DT FROM -X TO INFINITY 
EXNEI(X)=EXP(-X)*EI(X) (X .GT. 0) 

INTEGRAL(EXP(-T)/T) DT FROM X TO INFINITY* OFTEN 
DENOTED BY -EI (-X) =E1 (X) . (SEE AUTOMATIC COMPUTING 
METHODS FOR SPECIAL FUNCTIONS* PART II. THE EXPO- 
NENTIAL INTEGRAL EN(X)p J. OF RESEARCH NBSr 78Br 
OCTOBER-DECEMBER 1974* PP. 199-?16.) 
SHI(X) =INTEGRAL(SINH T/T)DT FROM TO X 
SHI(-X)=-SHI(X) 

CHKX)=GAMMA+LN X+INTEGRAL< (COSH T-n/T)DT FROM TO X 
CHI(-X)=CHI(X)-I PI 

GAMMA ( EULER • S CONSTANT ) = . 5772 1 56649 . . . 
SPECIAL CASES 
X=0 

SI(0)=SHI(0)=0 
CI(0)=EI(0)=EXNEI(O)=CHi(n)=-lNFINITY 

=-MAX. MACH. VALUE (RINF) 
LIMITING VALUES - X APPROACHES INFINITY 
SI(X)=PI/2 
CI(X)=0 

EI(X)=SHI(X)=CHI(X)=INFINTTY (RINF) 
EXNEI(X)=0 
USAGE. CALL SICIEI ( IC*X* SI * CI f CII *El rEXNEI 'SHI tCHl rCHII * 

IERR) 



FORMAL 
IC 



PARAMETERS 

INTEGER 
IC 
1 
2 
3 
4 
OR 



X REAL 

SI=SI(X) 

CI + I CII=CKX) 

EI=EI(X) 

EXNEI=EXP(-X)*EKX) 

SHI=SHI(X) 

CHI + I CHII=CHI(X) 

IERR INTEGER TYPE 

IERR=0 X .GE. 

IERR=1 X .LT. 



TYPE 

FUNCTIONS TO BE COMPUTED 
SI*CI 
ElfEXNEI 

ElrEX^EIfSHlfCHI 
SItCl'EI*EXNEI*SHI*CHI 
DOUBLE PRECISION TYPE 
(SAME TYPE AS X) 
t t 



INPUT 



t t 
t t 
t t 



Of 
Or 



INPUT 

OUTPUT 

OUTPUT 

OUTPUT 

OUTPUT 

OUTPUT 

OUTPUT 

OUTPUT 



NORMAL RETURN 
ERROR RETURN IF 
IC=2 



MODIFICATIONS. 
THE CODE IS 
WITH DOUBLE 

DOUBLE 
LARLYrFOR 

RINF 



SET UP FOR DOUBLE PRECISION COMPUTATION 

PRECISION TYPE STATEMENTS 

PRECISION FUNCTION REFERENCES AND*PARTICU- 
THE UNIVAC 1108 WITH (SEE DEFINITIONS BELOW) 
APPROX. 2**1023*ULSC=2**56*NBM=60 AND OTHER 
CONSTANTS IN DOUBLE PRECISION FORMAT TO 19 SIGNIFICANT 
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59* 

60* 

61* 

62* 

63* 

64* 

65* 

66* 

67* 

68* 

69* 

70* 

71* 

72* 

73* 

74* 

75* 

76* 

77* 

78* 

79* 

80* 

81* 

82* 

83* 

84* 

85* 

86* 

87* 

88* 

89* 

90* 

91* 

92* 

93* 

94* 

95* 

96* 

97* 

98* 

99* 

100* 

101* 

102* 

103* 

104* 

105* 

106* 

107* 

108* 

109* 

110* 

111* 

112* 

113* 

114* 

115* 

116* 



C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 



X=SUM OF SERIFS 
READILY BE ADAPTED TO 
FOR COMPLEX ARGUMENTS. 



FIGURES. ALL ABOVE ITEMS MUST BE CHANGED FOR SINGLE 
PRECISION COMPUTATIONS WITH DATA ADJUSTMENTS FOR OTHER 
COMPUTERS. 
AUXILIARY FUNCTIONS 

VARIOUS FUNCTIONS ARE AVAILABLE TO GREATER ACCURACY 
AT INTERMEDIATE POINTS IN THE SUBROUTINE* NAMELY* 

SI-(PI/2)=IMAG. PART OF THE CONTINUED FRACTION 
CI (EI AND CHD-GAMMA-LN 
CAUTION - THE SUBROUTINE CANNOT 
COMPUTE THE FUNCTIONS 
METHOD. T=ABS(X) 

POWER SERIES T .LE. PSLSC(=2) FOR SIfCI 

T .LE. AELL(=-LN(TOLER)) / FOR ElrSHIrCHI 
S I =SUMS ( SGN ( RK ) *TM ( RK ) ) 
CI=SUMC(SGN(RK)*TM(RK) ) 

+EULER+XLOG 
SHI=SUMOT(TM(RK)) 
CHI=SUMET(TM(RK)) 

+EULER+XLOG 
EI=SUMOT+SUMET+EULER+XLOG 
SGN(1)=1 

SGN(RK+1)=-SGN(RK) 
SGN(RK+1)=+SGN(RK) 
TM(RK)=( (T**RK)/(l*2t 
=PTM(RK)/RK 
PTM(1)=T 

PTM(RK+1)=PTM(RK)*(T/(RK+1) ) 
IF TM(RK)/SUM .LT. TOLER 



IP=-1 


RK=1*3* ...rRKO 


IP=+1 


RK=2#4r ... rRKE 


IP=-1 


RK=1*3, ...tRKO 


IP=+1 


RK=2f4r ...rRKE 



(X .GT. 0) 



RK=lf3r 
RK=2r4f 
•RK))/RK 



RK .GE. 1 



RKE= 


:RK 


WHERE 


SUM: 


=ABS(SUMC) 






IC = 


:1 OR 


4 








SUm: 


=SUMET 






IC= 


:2 OR 


3 










IC= 


:a. 


X 


.GT. 


PSLSC 


RKO: 


:RK 


WHERE 


SUM: 


=ABS(SUMS) 






IC= 


:1 OR 


4 








SUm: 


=SUMOT 






IC= 


:2 OR 


3 










IC= 


:'i« 


X 


.GT. 


PSLSC 



EXNEI= EI/EXP(T/2)/EXP(T/2) 
=(EI/EXPHT)/EXPHT 
CONTINUED FRACTION T .GT. PSLSC 
-CI+I(SI-PI/2>=E1(IT) 

=EXP(-IT)*(1 I/I 
1**2 I/I 
2**2 I/I 



(1+IT)- 
(3+IT)- 
(5+IT)-., 



) 



rRMF 



=EXP(-TT)*IKAM(RM) I/I BM(RM)) 

RM=! t2t . ., 
AM(1)=1 

AM(RM)r-(RM-l)**2 RM *GT. 1 
RM(RM)=2*RM-1. + IT=PMR+T RMI 
=EXP(-IT)*CFM/GM) 

=EXP(-IT)*(FMR+I H4I)/(GMR+I GMI ) 
=EXP(-IT>*F(RM) 
=(COST-I SINT)*(FP+I FT) 
-CI-H(SI-PI/2)=(FR*C0ST+FI*SINT)+ 

I(FI*COST-FR*SINT) 
IF RESQCRM) .LE. TOLSO (=T0LER**2) 
OR RFSO(RM) .GE. RESO(RM-l) 
(PESO .GE. RESOP) 
RMFrPM WHERE 

PES^=(MOn(l-F(RM-l)/F(RM)) ) **2 
ASYMPTOTIC EXPANSION T .GT. AFLL 
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117* 
118* 
119* 
120* 
121* 
122* 
123* 
12*4* 
125* 
126* 
127* 
128* 
129* 
130* 
131* 
132* 
133* 
134* 
135* 
136* 
137* 
138* 
139* 
140* 
141* 
142* 
143* 
144* 
145* 
146* 
147* 
148* 
149* 
150* 
151* 
152* 
153* 
154* 
155* 
156* 
157* 
158* 
159* 
160* 
161* 
162* 
163* 
164* 
165* 
166* 
167* 
168* 
169* 
170* 
171* 
172* 
173* 
174* 



C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 

c 
c 
c 
c 
c 
c 
c 
c 
c 



c 
c 



EI=(EXNEI*EXPHT)*EXPHT 

EXNEI= ( 1+SUME < TM < RK > ) ) /T RK=1 f 2 f . • . * RKF 

SHI=CHl=EI/2 

TM(RK)=(1*2...RK)/(T**RK) 
TM(0)=1 

TM(RK)=(RK/T)*TM(RK-1) RK .GE. 1 

IF TM(RK) .LI, TOLER (CONVERGENCE) RKF=RK OR 

TM(RK) .GE. TM(RK-l) (DIVERGENCE) RKF=RK-1 



RANGE. 
FOR 



SI(X)tCI(X)» ABS(X) .LT t 



FOR 
FOR 



FOR 



FOR 
ROUTINE) 



NRC=8> 
NBC=11) 
BIASED 
NUMBER 



ULSC (UPPER LIMIT 
SIN 9 COS 
X=APPR0XIMATELY 2**21 t MBM=27 
2**56 p MRM=60 
EXP(-X)*EKX) f X .LE. RTNF 

EI(X), X .LT. XMAXEI (APPROXIMATELY 92. 5* 

715. 6* 
NBC=NUMRER OF BINARY DIGITS IN THE 
CHARACTERISTIC OF A FLOATING POINT 
SHI(X) rCHl'(X) f ABS(X) .LT, XMAXHF 
X=APPR0XIMATELY ^3.2f NBC=8 
716.3* NBC=11 
ACCURACY. THE MAXIMUM RELATIVE ERROR* EXCEPT FOR REGIONS 
IN THE IMMEDIATE NEIGHBORHOOD OF ZFROS>ON THE 
UNIVAC 1108 IS 4.5(-7) FOR SINGLE PRECISION COM- 
PUTATION AND 7.5(-17) FOR DOUBLE PRECISION COM- 
PUTATION. 
PRECISION. VARIABLE - BY SETTING THE DESIRED VALUE OF NBM 

OR A PREDFTFRMINED VALUE OF TOLER 
MAXIMUM UNIVAC 1108 TIME/SHARING EXECUTIVE SYSTEM 
TIMING. NBM=27 NBM=60 
(SECONDS) .0093 .070 
STORAGE. 954 WORDS REQUIRED BY THE UNIVAC 1108 COMPILER 



SUBROUTINE SICIEI (ICf X*SI*CIf CI I >EI 9 EXNET rSHl ,CHI f 
1 CHII'IERR) 

MACHINE DEPENDENT STATEMENTS 
TYPE STATEMENTS 
DOUBLE PRECISION X f SI > CI * ClT 'EI rFXNEI ,SHT 9 CHI rCHI I 
DOUBLE PRECISION A t AELLf AM» AMINr ASUMSCt 

1 RMIfRMRrCOSTfEXPLfEXPHTf 

2 FI 9 F IP t FMI 9 FMM 1 1 * FMM1 R 9 FMM2 1 9 FMM2R * FMR p FR 9 FRP r 

3 GMIfGMMlItGMMlRfGMM2ItGMM2R»GMRf 

4 PSLL 9 PSLSC 9 PTM 9 RE 9 RESO # RESQP 9 RK 9 RM , 

5 SCC 9 SFM 1 9 SFMR f SGM 1 , SGMR 9 SGN 9 

6 SI NT f SUMC 9 SUME 9 SUMEO f SUMET r SUMOT 9 SUMS 9 SUMSC 9 

7 T 9 TEMP 9 TEMP A 9 TEMPR , TM 9 TM A X 9 TMM 1 r TOLER 9 TOLSQ 9 

8 XLOGtXMAXEIf XMAXHF 

DOUBLE PRECISION RlMFr ULSC »FULER *HALFPI 9 PI , AL0G2* 
1 ZERO fONF* TWO 9 FOUR 

DIMENSION A (4) 
EQUIVALENCE (FMRrA(l))r 



(FMTt A(2)) 9 (GMRf A(3) )r 



L (GMItA(4)) 

CONSTANTS 
DATA EULER/. 577215664901532860600/ 
DATA HALFPI/1.5707963267948Q6619D0/ 
DATA PI/3,14159265358979323^D0/ 
DATA ALOG2/.693147180559Q453094D0/ 
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175* DATA ZERO 9 ONE 9 TWO 9 FOUR / 

176* 1 0«0D0»l«DQ'2«D0'<t-*QG/ 

177* C RINF=MAXIMUM MACHINE VALUE 

178* C ULSC=MAXIMUM ARGUMENT FOR SIM , COS ROUTINE 

179* C APPROX. ?**(NRM-6) OR lO**(S-2) 

180* C (^SIGNIFICANT FIGURES) 

181* C NBM=ACCURACY DESIRED OR THE 

182* C MAXIMUM NUMBER OF BINARY DIGITS IN THE 

183* C MANTISSA OF A FLOATING POINT NUMBER 

184* C TOLER=UPPER LIMIT FOR RELATIVE ERRORS 

185* C r2**(~NRN)=:APPROX. 10** (-S) 

186* C TOLER PRECO^ p UTED MAY BE INSERTED IN A DATA STATEMENT AND 

187* C THE MBM DATA STATEMENT ELIMINATED 

188* DATA RlNF/. 8988465674311579^380308 / 

189* DATA ULSC/. 72057594037927936017/ 

190* DATA N8M / 60 / 

191* TOLER=TWO**(-NBM) 

192* C NOTE - ARGUMENT CHECKS PRECEDING FUNCTION REFERENCES 

193* C NECESSITATE ADDITIONAL MACHINE DEPENDENT STATEMENTS 

194* C IN THE STATEMENT NUMBER RANGE 140-150 

195* C INITIALIZATION OF OUTPUT FUNCTIONS 

196* SI=RINF 

197* CI=RINF 

198* CII=RINF 

199* EI=ZER0 

200* EXNEI=RINF 

201* SHI-ZERO 

202* CHI=ZER0 

203* CHII=RINF 

204* C VALIDITY CHECK ON INPUT PARAMETERS 

205* C INDICATOR CHECK 

206* C SET IND=IC 

207* C CHANGE IND=4 TF IC .LT. 1 OR .GT. 4 

208* IND=IC 

209* IF (IND .LT. 1) GO TO 10 

210* IF (IND .GT. 4) GO TO 10 

211* GO TO 20 

212* 10 IND=4 

213* C ARGUMENT CHECK 

214* C X .GE. TERP=0 

215* C X .LT. IERR=1 

216* C (ERPOR RETURN IF IC=?) 



217* 


20 


IERR=0 


218* 




T=X 


219* 


30 


IF (T) 40f50»90 


220* 


40 


T=-T 


221* 




IF (IND .EQ. 1) GO TO 30 


222* 




IERR=1 


223* 




IF (IND .NE. 2) GO TO 30 


224* 




IF (X .LT. ZERO) RETURN 


225* 


C 


SPECIAL CASES 


226* 


C 


X=0 


227* 


5n 


IF (IND-2) 80*70*60 


228* 


60 


SHI=ZER0 


229* 




CHI=-RINF 


230* 




CHII=ZER0 


231* 


70 


EI=-RINF 


232* 




EXNEI=-RINF 
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233* IF (IND .ME. 4 ) RETURN 

234* 80 5I=ZER0 

235* CI=-RINF 

236* CII=ZERO 

237* RETURN 

238* 90 IF (T # LT. ULSC) GO TO 140 

239* C ARS(X) .GF. ULSC 

240* IF (IND-2) 130*110*100 

241* 100 SHI=RINF 

242* CHI=RINF 

243* CHII=ZERO 

244* IF (IERR .EQ. 1) GO TO 120 

245* 110 EI=RINF 

246* EXNEI= ( ONE+ ( ONE/T ) ) /T 

247* 120 IF (INO .NE. 4) GO TO 1000 

248* 130 SI=HALFPI 

249* CI=ZERO 

250* CII=ZERO 

251* GO TO 1000 

252* C EVALUATIONS FOR ABS(X><=T) •GT. AMD .LT. ULSC 

253* C ADDITIONAL MACHINE DEPENDENT STATEMENTS 

254* C FUNCTION REFERENCES 

255* C CONTROL VARIABLES 

256* 140 XLOG=DLOG(T) 

257* 5INT=DSINCT) 

258* COST=DCOS(T) 

259* EXPL =DL0G(RINF) 

260* XMAXEI=EXPL+OLOG(EXPL+DLOG(FXPL) ) -ONE/EXPL 

261* XMAXHF=XMAXEI+AL0G2 

262* AELL=-OLOG ( TOLER ) 

263* AMIN=0NE/RINF 

264* PSLL=TWO*DSQRT ( AMIN ) 

265* PSLSC=TWO 

266* C EXPONENTIAL FUNCTION DETERMINATION 

267* IF (T .LE. TQLER) GO TO 150 

268* IF (T .GE. XMAXHF) 00 TO 160 

269* EXPHT=DEXP(T/TWO) 

270* GO TO 170 

27.1* 150 EXPHT=ONE 

272* GO TO 170 

273* 160 EXPHT=RINF 

274* C METHOD SELECTION 

275* 170 IF (T ,LE, PSLSC) GO TO 200 

276* IF (IND .EQ. 1) GO TO 500 

277* IF (IND .EQ. 4) GO TO 500 

278* 180 IF (T ,6T. AELL) GO TO 80Q 

279* GO TO 230 



280* 


C 


INDICATOR TO COMPUTE EI*SHI,CHI 


281* 


190 


IF (IND .EG. 1) GO TO 1000 


282* 




IND=3 


283* 




GO TO 180 


284* 


C 


METHOD Power SERIES 


285* 


C 


SI(X)fCKX) t T .LE, PSLSC 


286* 


C 


EKX)fSHI(y) rCHI(X)* T ,LE. AELL 


287* 


C 


LIMITING VALUES* T NEAR ZERO 


288* 


200 


IF (T .GT. PSLD GO TO 210 


289* 




SUMC=ZERO 


290* 




SUMET=ZERO 
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291* 




292* 




293* 




29'+* 


C 


295* 




296* 




297* 




298* 




299* 




300* 




301* 


C 


302* 




303* 




304* 




305* 




306* 


c 


307* 


c 


308* 




309* 




310* 




311* 


c 


312* 




313* 


c 


314* 




315* 




316* 


c 


317* 


c 


318* 




319* 




320* 




321* 




322* 




323* 




324* 




325* 




326* 


c 


327* 




328* 




329* 


c 


330* 




331* 


c 


332* 


c 


333* 




334* 




335* 




336* 




337* 




338* 




339* 




340* 




341* 




342* 




343* 




344* 


c 


345* 




346* 


c 


347* 




348* 


c 



INITIALIZATION FOR ST,CT 
GO TO 230 



SUMS=T 

SUMOTrT 
GO TO 360 

210 IF (IND .NEf 
220 SUMS=ZERO 
SUMC=ZERO 
SUMSC=ZERO 
S6N=0NE 
GO TO P40 

INITIALIZATION ^OR SHItCHHAND EI) 
230 SUMOT=ZER0 
SUMETrZERO 
SUMEO=ZERO 
IF (IND .EQ. 4) GO TO 220 

IP - INDICATOR FOR ODD OR 
EVEN TERMS 
240 IP=-1 
RK=ONE 
PTM=T 

COMPUTATION OF ( T**K ) / ( 1 *2. . .K ) /K 
250 TM=PTM/RK 

SUMMATION FOR SI (CI) 
IF (INO .NE. 1) GO TO 310 
260 SUMSC=SGN*TM+SUMSC 

RELATIVE ERROR FOR SI (CI) 
C PARTIAL SUM OF ALTERNATING OPD^VEN) TERMS MAY EOUAL ZERO 
ASUMSC=SUMSC 
270 IF (ASUMSC) 280'300r290 
280 ASUMSC=-ASUMSC 

GO TO 270 
290 RE=TM/ASUMSC 

GO TO 320 
300 RE=RINF 

GO TO 320 

SUMMATION Fop SHI (CHI) (AND ED 
310 SUMFO^TM+SUMEO 

IF (IMP .EQ. 4) GO TO 260 

RELATIVE ER p oR FOR SHI (CHI) 
RF=TM/SUMEO 

SIGN CHANGE AND SELECTION 
OF SUMS OF onD(FVEN) TERMS 
320 IF (IP .EQ. 1) GO TO 330 
SGN=-SGN 
SUMS-SUMSC 
SUMSC=SUMC 
SUMOT=SUMEO 
SUMEO=SUMET 
GO TO 340 
330 SUMC-SUMSC 
SUMSC-SUMS 
SUMET-SUMEO 
SUMEO^SUMOT 

RELATIVE ERROR CHECK 
340 IF (RF .LT. TOLER) GO TO ^60 

ADDITIONAL TERMS 
RKrRK+ONE 

UNDERFLOW Tf^ST 
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349* 
350* 
351* 
352* 
353* 
354* 
355* 
356* 
357* 
358* 
359* 
360* 
361* 
362* 
363* 
364* 
365* 
366* 
367* 
368* 
369* 
370* 
371* 
372* 
373* 
374* 
375* 
376* 
377* 
378* 
379* 
380* 
381* 
382* 
383* 
384* 
385* 
386* 
387* 
388* 
389* 
390* 
391* 
392* 
393* 
394* 
395* 
396* 
397* 
398* 
399* 
400* 
401* 
402* 
403* 
404* 
405* 
406* 



UNDERFLOWS AFFECTING ACCURACY ARE AVOIDED, ALL OTHFR 
UNDERFLOWS ARE ASSUMED TO BE SET EQUAL TO ZERO 
IF (T .GT. PSLSC) GO TO 350 
IF (PTM .LE. (AMIM*RK*RK)/T ) GO TO 360 
PTM=(T/RK)*PTM 
IP=-IP 
GO TO 250 

Sir CI EVALUATION 
IF (IND .NE. 1) GO TO 380 
SI=5UMS 

CI=(SUMC+XLOG)+EULER 
CII=ZERO 
GO TO 1000 

EI EVALUATION 
IF (X .LE. ZERO) GO TO 390 
EI= (SUMET-f SUMOT-f XLOG) +EULER 
EXNEI=(EI/EXPHT)/EXPHT 
IF (IND .EQ. 2) RETURN 

SHI t CHI EVALUATION 
SHI=SUM0T 

CHl=(EULER+SUMET)+XLOG 
CHII=ZER0 

4) GO TO 1000 



350 



360 
370 



380 



390 



C 
C 
C 
C 



500 



C 
C 
C 



510 



C 

C 

c 
c 
c 
c 



IF 
GO 



(IND .NE« 
TO 370 



METHOD CONTINUED FRACTION 

SI(X)fCI(X)f T .GT. PSLSC 
-CKT) + I (SI(T)-HALFPI)=F1(IT> 
INITIALIZATION 



SCC=RINF/F0UR 

TOLSQ=TOLFR*TOLER 

RM=ONE 

AM=ONE 

BMR=ONE 
BMI=T 

FMM2R=0NE 

FMM2I=ZERO 

GMM2R=ZERO 

GMM2I=ZERO 

FMM1R=ZER0 

FMM1I=ZER0 

GMM1R=0NE 

GMM1I=ZER0 

RESQP=RINF 

FRP=ZER0 

FIP=ZERO 



RECURRENCE RELATION 



FM=BM*FMM1 + AM*FMM2 
GM=BM*GMM1 + AM*GMM2 
FMR=PMR*FMM1R-BMI*FMM1I + AM*P : MM2R 
FMI=BMI*FMM1R+BMR*FMM1I+AM*FMM2I 
GMR=BMR*GMM1R-BMI*GMM1T+AM*GMM2R 
GMI=RMI*GMM1R+BMR*GMM1I+AM*GMM2I 

CONVFRGFNT F=FM/GM 

TESTS TO AVOID INCORRECT RESULTS 
DUE TO OVERFLOWS (UNDERFLOWS) 
FINDING MAXIMUM(=TMAX) OF 

ABSOLUTE OF FMR, GMR,FMI rGMl 
FOR SCALING PURPOSES 
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407* 
408* 
409* 
410* 
411* 
412* 
413* 
414* 
415* 
416* 
417* 
418* 
419* 
420* 
421* 
422* 
423* 
424* 
425* 
426* 
427* 
42R* 
429* 
430* 
431* 
432* 
433* 
434* 
435* 
436* 
437* 
438* 
439* 
440* 
441* 
442* 
443* 
444* 
445* 
446* 
447* 
448* 
449* 
450* 
451* 
452* 
453* 
454* 
455* 
456* 
457* 
458* 
459* 
460* 
461* 
462* 
463* 
464* 





TMAX=7ERO 




1 = 1 


520 


TEMPrA(T) 


530 


IF (TEMP) 540*560,550 


540 


TEMPs-TEMP 




GO TO 530 


550 


IF (TEMP .LE. TMAX) GO TO 560 




TMAX=TEMP 


560 


IF (I .GE. 4) GO TO 570 



1=1+1 

GO TO 520 
570 SFMR=FMR/TMAX 

SFMI=FMI/TMAX 

SGMR=GMR/TMAX 

SGMI-GMI/TMAX 

TEMP=SGMR*SGMR ♦ SGMI*SGMI 

FR=(SFMR*SGMR+SFMI*S6MI)/TEMP 

FI=(SFMI*SGMR-5FMR*SGMI)/TEMP 
C RELATIVE ERROR CHECK 

TEMP=FR*FR+FI*FI 

TEMPA=(FRP*FR+FIP*FI)/TE V 'P 

TEMPR=(FIP*FR-FRP*FI)/TEMP 

TEMP=ONE-TEMPA 

RESQ =TEMP*TEMP+TEMPB*TEMPR 

IF (RESQ .LE. TOLSQ) Go TO 590 

IF (RESQ .GE. RESQP) GO TO 580 
C ADDITIONAL CONVERGENTS 

AM=-RM*RM 

RM=RM+ONE 

BMR=RMR+TWO 

FMM2R=FMM1R 

FMM2I=FMM1I 

GMM2R=GMM1R 

GMM2I=GMM1I 

FMM1R=FMR 

FMM1I=FMI 

GMM1R=GMR 

GMM1I=GMI 

FRP=FR 

FIP=FI 

RESQP=RESQ 
C SCALING 

C SCALING SHOULD NOT BE DELETED AS THE VALUES OF FMR,FMI AND 
C GMR>GMI MAY OVERFLOW FOR SMALL VALUES OF T 

IF (TMAX .LT. SCC/CBMR-AM ) ) GO TO 510 

FMM2R=FMM2R/TMAX 

FMM2I=FMM2I/TMAX 

GMM2R=GMM2R/TMAX 

GMM2I=GMM2I/TMAX 

FMM1R=FMM1R/TMAX 

FMM1I=FMM1I/TMAX 

GMM1R=GMM1R/TMAX 

GMM1I=GMM1I/TMAX 

GO TO 510 



DIVERGENCE OF RELATIVE FRROR 
ACCEPT PRIOR CONVERGENT 



580 



FR=FRP 
FI=FIP 
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465* 


C 


SltCI EVALUATION 


466* 


590 


SI=FI*C0ST-FR*SINT+HALFPI 


467* 




CI=-(FR*COST+FI*SINT) 


468* 




CII=ZERO 


469* 




GO TO 190 


470* 


C 


METHOD ASYMPTOTIC EXPANSION 


471* 


C 


EKX)fEXNEI(X) X *GT. AELL 


472* 


C 


SHI(T)=CHim=FI(T)/2 T .GT. AELL 


473* 


C 


INITIALIZATION 


474* 


800 


IF (IND *NE. 2) GO TO 880 


475* 


810 


SUME=ZERO 


476* 




RK=ZERO 


477* 




TM=ONE 


478* 


C 


ADDITIONAL TF!RMS 


479* 


820 


TMM1=TM 


480* 




RK=RK+ONE 


481* 




TM=(RK/T)*TM 


482* 


C 


TOLERANCE CHECK 


463* 




IF CTM .LT. TOLER) GO TO 840 


484* 




IF (TM .GE. TMM1) GO TO 830 


485* 




SUME-SUME+TM 


486* 




GO TO 820 


487* 


C 


DIVERGENT PATH 


488* 


830 


SUME=SUME-TMM1 


489* 


C 


FXNEI EVALUATION 


490* 


840 


IF'CX .LT* ZERO) GO TO 870 


491* 




EXNEI=(ONE+SUME)/T 


492* 


C 


EI EVALUATION - X .l_T. XMAXEI 


493* 




IF (T .GE. XMAXFI) GO TO 850 


494* 




EI=(EXMEI*EXPHT)*EXPHT 


495* 




GO TO 860 


496* 


C 


EI - LIMITING VALUE* X .GE. XMAXEI 


497* 


850 


EI=RINF 


498* 


C 


SHI^CHI EVALUATION - T .LT. XMAXHF 


499* 


860 


IF (IND .EQ. 2) RETURN 


500* 


870 


IF (T ,GE. XMAXHF) 60 TO 1000 


501* 




SHI=( ( ( ( ONE+SUME)/T)/TWO)*EXPHT)*EXPHT 


502* 




CHI=SHI 


503* 




CHII=ZER0 


504* 




GO TO 1000 


505* 


C 


SHIrCHl - LIMITING VALUE 


506* 


C 


T .GE. XMAXHF 


507* 


880 


IF ( T .LT. XMAXHF) GO TO 810 


508* 




5HI=RINP 


509* 




CHI=RINF 


510* 




CHII=ZER0 


511* 




IF ( X .GT. ZERO) GO TO 810 


512* 




GO TO 1010 


513* 


C 


ADJUSTMENTS FOR X .LT. 


514* 


1000 


IF (X .GT. ZERO) RETURN 


515* 


1010 


IF (IC .EQ. 3) GO TO 1020 


516* 




SI=-SI 


517* 




Cllr-PI 


518* 




IF (IC .EQ. 1) RETURN 


519* 


1020 


SHI=-SHI 


520* 




CHII=-PI 


521* 




RETURN 


522* 




END 
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Table 1 



SI(X) 



CI(X) 



.0 

.1-001 
.2-001 
.3-001 
.4-001 
.5-001 
.6-001 
.7-001 
.8-001 
.9-001 



.0 

.9999944444 

. 1999955556 

.2999850004 

.3999644461 

.4999305607 

.5998800129 

.6998094724 

.7997156101 

.8995950984 



6111108276 
0888852607 
0499380108 
5106467200 
6366745212 
5920656146 
5377692911 
6294499144 
0144401781 



6470528517 
8665206276 
0675616964 
4469605623 
5334525894 
8561513856 
0744599342 
3834290267 
3419097919 



85973-002 
10625-001 
14691-001 
17391-001 
59473-001 
83629-001 
97997-001 
56205-001 
40853-001 



.4027979520 
.3334907338 
.2929567223 
.2642060133 
.2419141543 
.2237094916 
.2083269121 
.1950112552 
.1832754260 



9823920722 
8599613460 
9811175640 
3009493497 
5519082432 
8693029859 
9543103133 
,8007321930 
4358445302 



4397561453 
8102204879 
0020167602 
5931669690 
1138357349 
1937022506 
4741544991 
1863050699 
6218391310 



54440 + 001 
74545+001 
13443+001 
80028+001 
19704+001 
78885+001 
04089+001 
50094+001 
18002+001 



.1 + 000 
.2+000 
.3+000 
.4+000 
.5 + 000 
.6+000 
.7 + 000 
.8 + 000 
.9 + 000 



.9994446110 
. 1995560885 
.2985040438 
.3964614647 
.4931074180 
.5881288096 
.6812222391 
.7720957854 
.8604707107 



8276950160 
2623382140 
0704316138 
5137288302 
4306668916 
0808006689 
1661131088 
8199656025 
4529293277 



5921185541 
0456944764 
6446229574 
0334263135 
1626707572 
9647904006 
9506811453 
3889712479 
4085411696 



90930-001 
16595 + 000 
64345+000 
17445+000 
76465+000 
83682+000 
94252+000 
89549 + 000 
29011 + 000 



.1727868386 
. 1042205595 
.6491729329 
.3788093464 
.1777840788 
.2227070695 
.1005147070 
.1982786159 
.2760678304 



6572966389 
6727819753 
7116174495 
2524433208 
0661290133 
9279762526 
0889783268 
5246717701 
6777286015 



9772515290 
6291690667 
6181105135 
4661757229 
5810271070 
6027439020 
9135892810 
5833582233 
0434115538 



65736+001 
58895+001 
48965+000 
44184+000 
56908+000 
62561-001 
28616+000 
84097 + 000 
80334+000 



o 
so 



.1 + 001 
.2+001 
.3+001 
.4+001 
.5 + 001 
.6+001 
.7+001 
.8+001 
.9+001 

.1 + 002 

.2+002 
.3+002 
.4+002 
.5+002 
.6+002 
.7+002 
.8+002 
.9+002 

.1 + 003 
.2+003 
.3+003 
.4+003 
.5+003 
.6+003 
.7+003 
.8+003 
.9+003 

.1 + 004 



.9460830703 
.1605412976 
. 1848652527 
.1758203138 
.1549931244 
.1424687551 
.1454596614 
.1574186821 
.1665040075 

. 1658347594 
.1548241701 
. 1566756540 
.1586985119 
.1551617072 
.1586745616 
.1561594849 
.1572330886 
. 1575663406 

.1562225466 
.4568382339 
.1570881088 
.1572114869 
.1572565882 
.1572461233 
.1571993932 
.1571355087 
.1570721487 

.1570233121 



6718301494 
8026948485 
9994682563 
9490530581 
9446741372 
2805065357 
2480935906 
7069420520 
8296024951 

2188740493 
0434398401 
0303511109 
3547845067 
4859358947 
2599474123 
1780061055 
9124873153 
6574562607 

8890562933 
3394698333 
2137495192 
2738117518 
2431687035 
9493979398 
2374915706 
6214727479 
6829785964 

9687712181 



1353313823 


17966+000 


.3374039229 


0096813466 


2646203889 


15077+000 


7672014819 


85889 + 001 


.4229808287 


7486499569 


8565153198 


25589+000 


9773025111 


19732+001 


.1196297860 


0800032762 


6472281176 


67785+000 


0555930335 


85016+001 


-.1409816978 


8693041163 


9144898694 


03593+000 


7440840073 


06390+001 


-.1900297496 


5664387861 


8458900116 


30081+000 


6903102791 


71420 + 001 


-.6805724389 


3247126204 


1683048406 


17428-001 


1476849383 


61604+001 


.7669527848 


2184518382 


9157630314 


36847-001 


8297145120 


66585+001 


.1224338825 


3200955729 


2295958268 


73503+000 


0665342789 


71085+001 


.5534753133 


3133607085 


6416484497 


93252-001 


3097187938 


96725+001 


-.4545643300 


4455372634 


5328299526 


27853-001 


6364334212 


95137+001 


.4441982084 


5353316539 


7687169925 


70578-001 


8373130900 


67982 + 001 


-.3303241728 


2071143779 


2264409630 


03714-001 


7566596201 


46420+001 


. 1902000789 


6208766961 


9812034202 


67962-001 


2798559485 


93775+001 


-.5628386324 


1163054401 


8589549846 


40996-002 


2644013231 


99104+001 


-.4813243377 


4432152888 


7234385245 


11961-002 


2298220467 


79853+001 


.1092198847 


3464977037 


8113544132 


70120-001 


5125172966 


74798+001 


-.1240250115 


5070958192 


0458829452 


17856-001 


3805334080 


46545 + 001 


.9986124071 


6431336804 


5713422489 


45291-002 


5234513880 


45027+001 


-.5148825142 


6104921444 


4355390534 


44979-002 


5878557542 


35465 + 001 


-.4378446093 


0278256791 


6569771749 


32552-002 


5231225344 


08620+001 


-.3332199918 


5921117799 


7045482583 


87851-002 


0132144796 


40848+001 


-.2123988830 


8463489343 


4793666984 


71878-002 


3434162096 


10243+001 


-.9320008144 


0429025451 


7261962053 


86680-003 


3169426317 


07478+001 


.7641202377 


3809449584 


1487520193 


67238-004 


3702809228 


14464+001 


.7788100127 


3975633981 


0132017166 


41226-003 


0846382718 


80577+001 


.1118158760 


9180513118 


2640905181 


04941-002 


0335292159 


62388+001 


.1108585782 


0500159240 


9398971297 


71935-002 


4796277803 


63344+001 


.8263155110 


9068228200 


1773882343 


20723-003 


tt/2 = 1.57079 


63267 94896 61923 


13216 91639 


75144 







Table 2 





X 






EI(X) 






EXP(- 


-X)*EI(X) 






.0 
.1-001 


— X 

-.4017929465 


4266693867 


7534341058 


83082+001 


— X 

-.3977950399 


2615576971 


3408360703 


19476+001 




.2-001 


-.3314706894 


4101539023 


9227816770 


05873+001 


-.3249071300 


3015878880 


1888549128 


04425+001 




.3-001 


-.2899115723 


9402794440 


9618312164 


59357+001 


-.2813433905 


5380940257 


4977395174 


78007+001 




.4-001 


-.2601256577 


5728261871 


6902308916 


05352+001 


-.2499259848 


2574874044 


8538470778 


63194+001 




.5-001 


-.2367884598 


5793745242 


5950757147 


49992+001 


-.2252401503 


9907626374 


1216434311 


86807+001 




.6-001 


-.2175282915 


5516236358 


6147208749 


96411 + 001 


-.2048604300 


3782595071 


0809271777 


33217+001 




.7-001 


-.2010800063 


5428875184 


2790326756 


89313+001 


-.1874857552 


3138763197 


4860337705 


30829 + 001 




.8-001 


-.1866884102 


7729969797 


4868911203 


84179+001 


-.1723351232 


0791016360 


9393503175 


42803+001 




.9-001 


-.1738663750 


3469511928 


7135014650 


31175+001 


-.1589019022 


1427078812 


0803355460 


85013+001 




.1+000 


-.1622812813 


9692766749 


6568299922 


74753+001 


-.1468381756 


5476302940 


1938294633 


01893+001 




.2+000 


-.8217605879 


0240031565 


3310869899 


97772 + 000 


-.6728006649 


8313731896 


8600992150 


20087+000 




.3+000 


-.3026685392 


6582588446 


8136901536 


93942 + 000 


-.2242223687 


1524378927 


5487443901 


40039+000 




.4+000 


.1047652186 


1932479322 


8929763883 


46026+000 


.7022622616 


7839611020 


0870681109 


89467-001 




.5+000 


.4542199048 


6317357992 


0523812662 


80237+000 


.2754982985 


5127026213 


2219977961 


61145+000 




.6+000 


.7698812899 


3735943709 


3137330257 


85125+000 


.4225198103 


2870176143 


1329955834 


52537+000 




.7 + 000 


.1064907194 


6242905405 


7162359274 


67593+001 


.5288172627 


5216097535 


3053374545 


58451 + 000 




.8+000 


.1347396548 


2123259381 


1899638330 


35829 + 001 


.6054242952 


6336443347 


4212274854 


28025+000 




.9 + 000 


.1622811713 


6968674413 


4419142440 


54222 + 001 


.6597860062 


6079393891 


3476848207 


67092+000 




.1 + 001 


.1895117816 


3559367554 


6652093433 


16343 + 001 


.6971748832 


3506606876 


5478681919 


55160+000 




.2+001 


.4954234356 


0018901633 


7950513022 


70353+001 


.6704827097 


9007328104 


3223808083 


90335+000 




.3+001 


.9933832570 


6254165580 


0833601921 


67653+001 


.4945764013 


4864123502 


8769696756 


17092+000 


OO 


.4+001 


. 1963087447 


0056220022 


6457202797 


23839+002 


.3595520078 


6362069617 


7421675826 


60745+000 


o 


.5 + 001 


.4018527535 


5803177455 


0914217937 


95867+002 


.2707662554 


9105719558 


0378348094 


65324+000 




.6+001 


.8598976214 


2439204803 


5834003079 


90690+002 


.2131473100 


8159360315 


3856832629 


71120+000 




.7+001 


.1915047433 


3550139595 


3063148272 


45695+003 


.1746297217 


6579015129 


0938352318 


25166+000 




.8+001 


.4403798995 


3483826899 


7424596659 


39339+003 


. 1477309983 


7340099664 


7174203305 


76383+000 




.9+001 


. 1037878290 


7170895876 


5757322679 


36222+004 


. 1280843565 


2321386813 


5996461057 


99933+000 




.1 + 002 


.2492228976 


2418777591 


3844014399 


85248+004 


.1131470204 


7341077803 


4051681354 


47701 + 000 




.2+002 


.2561565266 


4056588820 


4811208040 


98072+008 


.5279779527 


9648132254 


5779817983 


87276-001 




.3+002 


.3689732094 


0727419706 


4006328910 


84575+012 


.3452712179 


2361846131 


5818701748 


98032-001 




.4+002 


.6039718263 


6112415783 


5923141851 


06913+016 


.2565886278 


5975145205 


8133890489 


43081-001 




.5 + 002 


. 1058563689 


7131690963 


0615414332 


29987+021 


.2041704555 


5943987333 


5659940616 


10444-001 




.6 + 002 


.1936182213 


9292765388 


2072596687 


52373+025 


.1695420039 


4813286610 


0796814766 


65768-001 




.7+002 


.3646352759 


5797356367 


1345154856 


77168+029 


.1449589211 


5101025543 


4076487427 


26111-001 




.8+002 


.7014600004 


9047999696 


2996948047 


97919 + 033 


.1266031055 


4032884299 


9914561805 


60367-001 




.9+002 


.1371416869 


5072519995 


0222595544 


66918+038 


.1123740714 


1495976244 


6676886934 


83357-001 




.1 + 003 


.2715552744 


8538798219 


1401464231 


08254+042 


.1010206252 


7748357112 


3003599184 


99944-001 




.2+003 


.3631235233 


1593568523 


9671004384 


64250+085 


.5025253826 


9333012313 


0841542901 


12999-002 




.3+003 


.6496482508 


0886657890 


2569189493 


42473+128 


.3344519269 


3037826333 


1099042135 


67540-002 




.4+003 


. 1308647281 


7074277342 


4939487839 


11718+172 


.2506281486 


7484941759 


1294495095 


34843-002 




.5+003 


.2812821397 


8862943374 


7493151789 


64387+215 


.2004016096 


7757734712 


6302843056 


26131-002 




.6+003 


.6298882891 


3879314245 


2569543511 


94621 + 258 


. 1669453750 


3112400160 


6469063742 


31301-002 




.7 + 003 


.1450978736 


0525608526 


2088252210 


93011 + 302 


.1430618100 


9351634011 


0181288435 


03898-002 




.8+003 


.3412238865 


4483770461 


9667674034 


14661 + 345 


.1251566420 


9721409147 


4398485683 


.85823-002 




.9+003 


.8152195006 


2752522682 


3560124386 


87339+388 


.1112348431 


6823895903 


3663223603 


11253-002 




.1 + 004 


.1972045137 


1412383028 


0964504841 


20236+432 


.1001002006 


0241207250 


8068654920 


21171-002 























Table 3 



SHI(X) 



CHI(X) 



.1-001 


. 1000005555 


5722222505 


6692404330 


26380-001 


-.4027929520 


9823916092 


8101265102 


13346+001 


.2-001 


.2000044444 


9777814059 


1136870933 


39861-001 


-.3334707338 


8599317164 


5139185479 


39271 + 001 


.3-001 


.3000150004 


0500619903 


9860097325 


50634-001 


-.2929117223 


9807800640 


0016913137 


84864+001 


.4-001 


.4000355572 


6226866293 


4193925798 


75059-001 


-.2641260133 


2990530534 


6244248174 


04103 + 001 


.5-001 


.5000694496 


5299922650 


1974168087 


39980-001 


-.2417891543 


5446744469 


0970498828 


37392 + 001 


.6-001 


.6001200129 


6079350024 


5724232696 


91766-001 


-.2235294916 


8477029858 


8604451076 


93328+001 


.7-001 


.7001905835 


6955665134 


1973404093 


50990-001 


-.2080819121 


8998431835 


6210060797 


82823+001 


.8-001 


.8002844990 


6372249715 


1895071106 


72843-001 


-.1946912552 


6793692294 


6387861914 


90908+001 


.9-001 


.9004050984 


2855835469 


1051467784 


84646-001 


-.1828704260 


1898070283 


4045529328 


16021+001 


.1 + 000 


. 1000555722 


2505699555 


7615329532 


17784+000 


-.1722868386 


1943336705 


2329832875 


96531 + 001 


.2 + 000 


.2004449781 


4074638634 


0730853837 


22252+000 


-.1022205566 


0431467019 


9404172373 


72002+001 


.3+000 


.3015040562 


0501041398 


1095310302 


38247+000 


-.6041725954 


7083629844 


9232211839 


32189 + 000 


.4+000 


.4035726687 


4249363590 


5979378947 


55253+000 


-.2988074501 


2316884267 


7049615064 


09227 + 000 


.5 + 000 


.5069967498 


1966719583 


3659875988 


94380+000 


-.5277684495 


6493615913 


1360633261 


41435-001 


.6 + 000 


.6121303965 


6338077262 


4562784597 


54146+000 


. 1577508933 


7397866446 


8574545660 


30978+000 


.7 + 000 


.7193380189 


2889984241 


9121259127 


50575 + 000 


.3455691756 


9539069815 


2502333619 


25356+000 


.8+000 


.8289965633 


7893448638 


6910469189 


60092 + 000 


.5183999848 


3339145173 


2085914113 


98201 + 000 


.9+000 


.9414978265 


1143354092 


2701645733 


42970+000 


.6813138871 


8543390042 


1489778671 


99251 + 000 


.1 + 001 


. 1057250875 


3757285145 


7184235489 


58780+001 


.8378669409 


8020824089 


4678579435 


75631 + 000 


.2 + 001 


.2501567433 


3549756414 


7337248272 


75424+001 


.2452666922 


6469145219 


0613264749 


94929+001 


.3+001 


.4973440475 


8598067977 


1041838252 


27051 + 001 


.4960392094 


7656097602 


9791763669 


40601 + 001 


.4+001 


.9817326911 


2330344645 


6229756992 


81526+001 


.9813547558 


8231855580 


8342270979 


56862+001 


.5+001 


.2009321182 


5697226390 


4443761778 


82843+002 


.2009206353 


0105951064 


6470456159 


13024+002 


.6+001 


.4299506111 


2445683731 


1213478510 


53231 + 002 


.4299470102 


9993521072 


4620524569 


37459+002 


.7 + 001 


.9575242940 


8616503145 


6397896419 


56496+002 


.9575231392 


6884892807 


4233586305 


00452+002 


.8+001 


2201899686 


0023055646 


1163184608 


69467+003 


.2201899309 


3460771253 


6261412050 


69872+003 


.9 + 001 


.5189391515 


8222188283 


1922673971 


09373+003 


.5189391391 


3486770482 


5650552822 


52849+003 


.1 + 002 


.1246114490 


1994233444 


1188221070 


06923 + 004 


.1246114486 


0424544147 


2655793329 


78325+004 


.2+002 


. 1280782633 


2028294459 


4181868552 


98444+008 


.1280782633 


2028294361 


0629339487 


99627+008 


.3+002 


.1844866047 


0363709853 


2003165966 


19888+012 


.1844866047 


0363709853 


2003162944 


64687+012 


.4+002 


.3019859131 


8056207891 


7961570925 


53457+016 


.3019859131 


8056207891 


7961570925 


53456+016 


.5 + 002 


.5292818448 


5658454815 


3077071661 


49936+020 


.5292818448 


5658454815 


3077071661 


49936+020 


.6+002 


.9680911069 


6463826941 


0362983437 


61866 + 024 


.9680911069 


6463826941 


0362983437 


61866+024 


.7+002 


.1823176379 


7898678183 


5672577428 


38584+029 


.1823176379 


7898678183 


5672577428 


38584+029 


.8+002 


.3507300002 


4523999848 


1498474023 


98960+033 


.3507300002 


4523999848 


1498474023 


98960+033 


.9+002 


.6857084347 


5362599975 


1112977723 


34589+037 


.6857084347 


5362599975 


1112977723 


34589+037 


.1 + 003 


.1357776372 


4269399109 


5700732115 


54127+042 


.1357776372 


4269399109 


5700732115 


54127+042 


.2+003 


.1815617616 


5796784261 


9835502192 


32125 + 085 


.1815617616 


5796784261 


9835502192 


32125+085 


.3+003 


.3248241254 


0443328945 


1284594746 


71236+128 


.3248241254 


0443328945 


1284594746 


71236+128 


.4+003 


.6543236408 


5371386712 


4697439195 


58592+171 


.6543236408 


5371386712 


4697439195 


58592+171 


.5 + 003 


.1406410698 


9431471687 


3746575894 


82193+215 


.1406410698 


9431471687 


3746575894 


82193+215 


.6+003 


.3149441445 


6939657122 


6284771755 


97311 + 258 


.3149441445 


6939657122 


6284771755 


97311 + 258 


.7 + 003 


.7254893680 


2628042631 


0441261054 


65053 + 301 


.7254893680 


2628042631 


0441261054 


65053+301 


.8+003 


.1706119432 


7241885230 


9833837017 


07331 + 345 


.1706119432 


7241885230 


9833837017 


07331 + 345 


.9+003 


.4076097503 


1376261341 


1780062193 


43669+388 


.4076097503 


1376261341 


1780062193 


43669+388 


.1 + 004 


.9860225685 


7061915140 


4822524206 


01178+431 


.9860225685 


7061915140 


4822524206 


01178+431 





